#pragma once
#include <cstddef>
#include <cstdio>
#include <tuple>
#include <cctype>
#include <type_traits>
#include <cassert>
#include "cell.h"
#include "esmd_types.h"
#include "list_cpp.hpp"
#include "htab_cpp.hpp"
#include "utils.hpp"
#include "unitconv.hpp"
#include "tuple_hash.hpp"
#include "param_mapper.hpp"
#include "listed.hpp"
#include "treduce.hpp"
template <typename T, int IFunc>
struct gmxfunc {
  typedef T param_t;
  static const int ifunc = IFunc;
};
template <typename... Ts>
struct gmxfunclist {
};
template <>
struct gmxfunclist<> {
};
template <typename T, typename... Ts>
struct gmxfunclist<T, Ts...> {
  typedef typename T::param_t param_t;
  static constexpr int ifunc = T::ifunc;
};
template <typename... Ts>
struct gmxvparam {
  static constexpr size_t size = treduce::max(sizeof(typename Ts::param_t)...);
  static constexpr size_t align = treduce::max(alignof(typename Ts::param_t)...);
  typedef utils::type_list<typename Ts::param_t...> variants;
  static gmxfunclist<Ts...> funclist;

  alignas(align) unsigned char data[size];
  template <typename T>
  T &as() {
    static_assert(variants::template has<T>, "conversion type T not in type list Ts...");
    return *(T *)data;
  }

  template <typename T>
  operator T &() {
    return *(T *)(data);
  }
};

struct gmx_pair_route;
struct gmxtop {
  enum directive {
    DIR_DEFAULTS, //0
    DIR_ATOMTYPES, //1
    DIR_BONDTYPES, //2
    DIR_PAIRTYPES, //3
    DIR_ANGLETYPES,//4
    DIR_DIHEDRALTYPES,//5
    DIR_CONSTRAINTTYPES,//6
    DIR_NONBONDED_PARAMS,//7
    DIR_MOLECULETYPE,//8
    DIR_ATOMS,//9
    DIR_BONDS,//10
    DIR_PAIRS,//11
    DIR_ANGLES,//12
    DIR_DIHEDRALS,//13
    DIR_CONSTRAINTS,//14
    DIR_EXCLUSIONS,//15
    DIR_SYSTEM,//16
    DIR_MOLECULES,//17
    DIR_SETTLES,//18
    DIR_OTHERS,//19
    DIR_FILEEND,
    DIR_BROKEN
  };
  struct default_param {
    int nbfunc;
    int comb_rule;
    bool gen_pairs;
    real fudge_lj;
    real fugde_qq;
    DEFTIE(tie, nbfunc, comb_rule, gen_pairs, fudge_lj, fugde_qq)
  };
  struct atom_type {
    elemtag name;
    int atnum;
    real mass, charge;
    char ptype;
    real sigma, epsilon;
    DEFTIE(tie, name, atnum, mass, charge, ptype, sigma, epsilon);
  };

  typedef gmxvparam<gmxfunc<harmonic_bond_param, 1>,
                    gmxfunc<harmonic_bond_param, 2>>
      bond_param;

  struct bond_type {
    elemtag i, j;
    int ifunc;
    bond_param param;
    DEFTIE(tie, i, j, ifunc)
    bool operator<(const bond_type &o) const {
      return tie() < o.tie();
    }
  };

  typedef gmxvparam<gmxfunc<harmonic_angle_param, 1>,
                    gmxfunc<urey_bradly_param, 5>>
      angle_param;

  struct angle_type {
    elemtag i, j, k;
    int ifunc;
    angle_param param;
    DEFTIE(tie, i, j, k, ifunc)
    bool operator<(const angle_type &o) const {
      return tie() < o.tie();
    }
  };
  struct distance_constraint {
    real r0;
    DEFTIE(tie, r0)
  };
  typedef gmxvparam<gmxfunc<distance_constraint, 1>,
                    gmxfunc<distance_constraint, 2>>
  constraint_param;
  struct constraint_type {
    elemtag i, j;
    int ifunc;
    constraint_param param;
    DEFTIE(tie, i, j, ifunc)
    bool operator<(const constraint_type &o) const {
      return tie() < o.tie();
    }
  };
  typedef gmxvparam<gmxfunc<cosine_torsion_param, 1>,
                    gmxfunc<cosine_torsion_param, 9>,
                    gmxfunc<ryckaert_bellmans_param, 3>,
                    gmxfunc<harmonic_improper_param, 2>>
      dihedral_param;

  struct dihedral_type {
    elemtag i, j, k, l;
    int ifunc;
    dihedral_param param;
    DEFTIE(tie, i, j, k, l, ifunc)
    bool operator<(const dihedral_type &o) const {
      return tie() < o.tie();
    }
  };
  // struct lj_param {
  //   real sigma, epsilon;
  //   DEFTIE(tie, sigma, epsilon)
  // };
  typedef gmxvparam<gmxfunc<ljcoul_param, 1>> pair_param;
  typedef gmxvparam<gmxfunc<routed_ljcoul_param, 1>> routed_pair_param;

  struct pair_type {
    elemtag i, j;
    int ifunc;
    pair_param param;
    DEFTIE(tie, i, j, ifunc)
    bool operator<(const pair_type &o) const {
      return tie() < o.tie();
    }
  };
  struct bond {
    typedef bond_type param_type;
    tagint i, j;
    int ifunc;
    bool override_param;
    bond_param param;
    int pid;
    DEFTIE(tie_tags, i, j)
    static constexpr int nord = 2;
    std::tuple<tagint, tagint> tpl_tags(int ord) {
      if (ord == 0) return std::make_tuple(i, j);
      if (ord == 1) return std::make_tuple(j, i);
      return std::make_tuple(-1, -1);
    }
  };
  struct angle {
    typedef angle_type param_type;
    tagint i, j, k;
    int ifunc;
    bool override_param;
    angle_param param;
    int pid;
    DEFTIE(tie_tags, i, j, k)
    static constexpr int nord = 2;
    std::tuple<tagint, tagint, tagint> tpl_tags(int ord) {
      if (ord == 0) return std::make_tuple(i, j, k);
      if (ord == 1) return std::make_tuple(k, j, i);
      return std::make_tuple(-1, -1, -1);
    }
  };
  struct dihedral {
    typedef dihedral_type param_type;
    static constexpr int impr_ifuncs = 1<<4 | 1<<5;
    tagint i, j, k, l;
    int ifunc;
    bool override_param;
    dihedral_param param;
    int pid;
    DEFTIE(tie_tags, i, j, k, l)
    static constexpr int nord = 6;
    std::tuple<tagint, tagint, tagint, tagint> tpl_tags(int ord) {
      switch (ord) {
        case 0:          return std::make_tuple(i, j, k, l);
        case 1:          return std::make_tuple(l, k, j, i);
        case 2:          return std::make_tuple(-1, j, k, -1);
        case 3:          return std::make_tuple(-1, k, j, -1);
        case 4:          return std::make_tuple(i, -1, -1, l);
        case 5:          return std::make_tuple(l, -1, -1, i);
        default:
        return std::make_tuple(-1, -1, -1, -1);
      }
    }
  };
  struct pair {
    tagint i, j;
    int ifunc;
    bool override_param;
    pair_param param;
    int pid;
    DEFTIE(tie_tags, i, j);
  };
  struct routed_pair {
    tagint i, j, k;
    int ifunc;
    bool override_param;
    routed_pair_param param;
    typedef pair_type param_type;
    int pid;
    DEFTIE(tie_tags, i, j, k);
    static constexpr int nord = 2;
    std::tuple<tagint, tagint> tpl_tags(int ord) {
      switch (ord) {
      case 0:          return std::make_tuple(i, j);
      case 1:          return std::make_tuple(j, i);
      default:
        return std::make_tuple(-1, -1);
      }
    }
  };

  struct constraint {
    tagint i, j;
    int ifunc;
    bool override_param;
    constraint_param param;
    int pid;
    DEFTIE(tie_tags, i, j);
    static constexpr int nord = 2;
    std::tuple<tagint, tagint> tpl_tags(int ord) {
      switch (ord) {
        case 0:          return std::make_tuple(i, j);
        case 1:          return std::make_tuple(j, i);
        default:
          return std::make_tuple(-1, -1);
      }
    }
    bool operator<(const constraint &o){
      return i < o.i;
    }
  };
  struct atom {
    elemtag type, res_name, atom_name, atom_res_name;
    tagint atom_num, res_num, cg_num;
    real q, mass;
    //nr  type  resnr residue atom cgnr charge  mass
    DEFTIE(tie, atom_num, type, res_num, res_name, atom_res_name, cg_num, q, mass)
  };
  struct exclusion {
    tagint i, j;
    DEFTIE(tie, i, j);
    bool operator<(const exclusion &e) const {
      return tie() < e.tie();
    }
  };
  struct settle {
    int i, j;
    real r0oh, r0hh;
    bool present;
  };
  struct molecule_type {
    const char *name;
    int nrexcl;
    esmd::list<bond> bonds;
    esmd::list<angle> angles;
    esmd::list<dihedral> dihedrals;
    esmd::list<atom> atoms;
    esmd::list<pair> pairs;
    esmd::list<routed_pair> routed_pairs;
    esmd::list<exclusion> excls;
    esmd::list<constraint> constraints;
    settle settles;
    molecule_type() : bonds("moltype/bonds"), angles("moltype/angles"), dihedrals("moltype/dihedrals"), atoms("moltype/atoms"), pairs("moltype/pairs"), excls("moltype/excls"), constraints("moltype/constraints"), routed_pairs("moltype/routed pairs") {
      settles.present = false;
    }
    template<typename ...Ts, int ...Is>
    auto get_types(std::tuple<Ts...> ids, utils::iseq<int, Is...> seq) const {
      elemtag X = str2tag("X");
      return std::make_tuple(std::get<Is>(ids) == -1 ? X : atoms[std::get<Is>(ids) - 1].type...);
    }
    template<typename ...Ts>
    auto get_types(std::tuple<Ts...> ids) const {
      return get_types(ids, utils::make_iseq<int, sizeof...(Ts)>{});
    }
    void build_exclusions(gmxtop &top);
    template<typename T>
    esmd::list<T> &get_tops(){
      constexpr int ord = utils::get_ord<utils::type_list<bond, angle, dihedral, pair, routed_pair>, T>::value;
      return std::get<ord>(std::tie(bonds, angles, dihedrals, pairs, routed_pairs));
    }
    void find_clusters(listed_force<rigid_param> &rigids, param_mapper<rigid_param> &param_mapper, esmd::htab<tuple_key<elemtag, elemtag>, double> &r0s);
    void route_pairs();
    void sort_rigids();
    void convert_settle();
    //void export_special_pairs(listed_force<special_pair> &excls, listed_force<special_pair> &scals);
    //void fudge_pairtypes(esmd::list<atom_type> &atom_types, int comb_rule, double fudge_lj, double fudge_qq);
  };

  struct molecule {
    const char *name;
    tagint cnt;
    DEFTIE(tie, name, cnt)
  };
  esmd::htab<cchar_key, directive> dirtab;
  esmd::list<atom_type> atom_types;
  esmd::list<bond_type> bond_types;
  esmd::list<angle_type> angle_types;
  esmd::list<dihedral_type> dihedral_types;
  esmd::list<pair_type> pair_types;
  esmd::list<constraint_type> constraint_types;
  esmd::list<molecule_type> molecule_types;
  esmd::list<molecule> molecules;

  char *systemname;
  default_param defaults;
  gmxtop();
  void parse_file(FILE *f);

  directive parse_directive(const char *line);
  directive parse_others(FILE *f);
  directive parse_system(FILE *f);
  directive parse_moleculetype(FILE *f);
  template <typename T>
  directive parse_types(esmd::list<T> &list, FILE *f);
  template <typename T>
  directive parse_topologies(esmd::list<T> &list, FILE *f);
  template <typename T>
  directive parse_to_list(esmd::list<T> &list, FILE *f);
  directive parse_defaults(FILE *f);
  directive parse_settles(settle &set, FILE *f);
  directive parse_exclusions(esmd::list<exclusion> &list, FILE *f);
  
  // void make_mappers();
  void find_parameters();
  void dump(FILE *f);
  /*Convert GROMACS units to NAMD units: Angstrom, KCal/Mole, g/Mole*/
  void unitconv(unit lunit, unit eunit, unit runit);
  template<typename T>
  auto &get_params(){
    constexpr int ord = utils::get_ord<utils::type_list<bond_type, angle_type, dihedral_type, pair_type>, T>::value;
    return std::get<ord>(std::tie(bond_types, angle_types, dihedral_types, pair_types));
  }
  auto &get_moltype(const char *name){
    for (int i = 0; i < molecule_types.size(); i ++){
      if (!strcmp(name, molecule_types[i].name)) {
        return molecule_types[i];
      }
    }
    return *(molecule_type*)(NULL);
  }

  template<typename ParamAdapter, typename ToParam>
  void export_listed_adapt(listed_force<ToParam> &listed, param_mapper<ToParam> &mapper, int ifunc_sels);
  template<typename ToParam>
  void export_listed(listed_force<ToParam> &listed, param_mapper<ToParam> &mapper, int ifunc_sels);
  void export_excl_scal(cellgrid_t &grid);
  void export_special_pairs(listed_force<special_pair> &excls);
  void build_exclusions();
  void find_clusters(listed_force<rigid_param> &rigids, param_mapper<rigid_param> &param_mapper);
  void route_pairs();
  void fudge_pairtypes();
  void create_vdw14_params(esmd::htab<tag_key, std::tuple<real, real, real, real>> &map);
  void sort_rigids();
  void convert_settle();
  bool verbose = false;
};

template <typename ParamType>
using gmx_top_of = typename utils::get_type<
  utils::type_list<
    gmxtop::bond,
    gmxtop::angle,
    gmxtop::angle,
    gmxtop::dihedral,
    gmxtop::dihedral,
    gmxtop::dihedral,
    gmxtop::dihedral,
    gmxtop::routed_pair>,
  utils::get_ord<utils::type_list<
                   harmonic_bond_param,
                   harmonic_angle_param,
                   urey_bradly_param,
                   ryckaert_bellmans_param,
                   cosine_torsion_param,
                   extended_ryckaert_bellmans_param,
                   harmonic_improper_param,
                   routed_ljcoul_param>,
                 ParamType>::value
  >::type;
